function current = ghk(voltage, out, in, z)

if nargin < 4, z = 2; end
if nargin < 3, in = 75e-9; end
if nargin < 2, out = 1.5e-3; end

RTF = 8.314 * (37 + 273.15) / 96485e-3; % in mV
useVoltage = voltage;
useVoltage(useVoltage==0) = 1e-6;

current = useVoltage .* (in - out*exp(-z*useVoltage/RTF)) ./ (1 - exp(-z*useVoltage/RTF));
